General Recipe for Designing Photonic 
Crystal Cavities 



Dirk Englund 

Department of Applied Physics, Stanford University, Stanford, CA 94305 

Ilya Fushman 

Department of Biophysics, Stanford University, Stanford, CA 94305 

Jelena Vuckovic 

Department of Electrical Engineering, Stanford University, Stanford, CA 94305 
All authors contributed equally. 

Abstract: We describe a general recipe for designing high-quality 
factor (Q) photonic crystal cavities with small mode volumes. We first 
derive a simple expression for out-of-plane losses in terms of the £-space 
distribution of the cavity mode. Using this, we select a field that will result 
in a high Q. We then derive an analytical relation between the cavity field 
and the dielectric constant along a high symmetry direction, and use it to 
confine our desired mode. By employing this inverse problem approach, 
we are able to design photonic crystal cavities with Q > 4 • 10 6 and mode 
volumes V ~ (X /«) 3 . Our approach completely eliminates parameter space 
searches in photonic crystal cavity design, and allows rapid optimization of 
a large range of photonic crystal cavities. Finally, we study the limit of the 
out-of-plane cavity Q and mode volume ratio. 
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1. Introduction 

One of the most interesting applications of photonic crystals (PhCs) is the localization of light to 
very small mode volumes - below a cubic optical wavelength. The principal confinement mech- 
anism is Distributed Bragg Reflection (DBR), in contrast to, e.g., microspheres or microdisks, 
which rely solely on Total Internal Reflection (TIR). However, since three-dimensional (3D) 
PhCs (employing DBR confinement in all 3D in space) have not been perfected yet, the main- 
stream of the PhC research has addressed 2D PhCs of finite depth, employing DBR in 2D and 
TIR in the remaining ID. Such structures are also more compatible with the present microfab- 
rication and planar integration techniques. These cavities can still have small mode volumes, 
but the absence of full 3D confinement by DBR makes the problem of the high-quality fac- 
tor (high-Q) cavity design much more challenging. The main problem is out-of-plane loss by 
imperfect TIR, which becomes particularly severe in the smallest volume cavities. 

Over the past few years, many approaches have been proposed to address this issue, but they 
all focused on the optimization of a particular cavity geometry and a particular mode supported 
by it [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Some of the proposed cavities seem promising and have 
already been proven useful as components of lasers or add/drop filters. However, a general 
recipe for designing optimized nanocavities is missing, and it is also not known if cavities 
better than the ones presently known are possible; these are exactly the issues that we will 
address in this article. In Section 2, we will first derive a simple set of equations for calculation 
of cavity Q and mode volume. In Section 3 we estimate the optimum £-space distribution of 
the cavity mode field. In Section 4, we will finally address the question of finding PhC cavities 
with maximum possible figures of merit for various applications. In this process, we start from 
the optimum fc-space distribution of the cavity field; then we derive an approximate analytical 
relation between the cavity mode and the dielectric constant along a direction of high symmetry, 
and use it to create a cavity that supports the selected high-Q mode in a single step. Thus, we 
eliminate the need for trial and error or other parameter search processes, that are typically used 
in PhC cavity designs. Furthermore, we study the limit of out-of-plane Q factor for a given V 
of the cavity mode with a particular field pattern. 

2. Simplified relation between Q of a cavity mode and its fc-space Distribution 

In order to simplify PhC cavity optimization, we derive an analytical relation between the near- 
field pattern of the cavity mode and its quality factor in this section. Q measures how well the 
cavity confines light and is defined as 



(1) 



where co is the angular frequency of the confined mode. The mode energy is 



(U) = J l -(eE 2 + pLH 2 )dV (2) 

The difficulty lies in calculating P, the far-field radiation intensity. 

Following our prior work [2], we consider the in-plane and out-of-plane mode loss mecha- 
nisms in two-dimensional photonic crystals of finite depth separately: 

(P) = (P n } + (P x ) (3) 

(4) 
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In-plane confinement occurs through DBR. For frequencies well inside the photonic band gap, 
this confinement can be made arbitrarily high (i.e., Qn arbitrarily large) by addition of PhC 
layers. On the other hand, out-of-plane confinement, which dictates Q±, depends on the modal 
k-distribution that is not confined by TIR. This distribution is highly sensitive to the exact mode 
pattern and must be optimized by careful tuning of the PhC defect. Assuming that the cavity 
mode is well inside the photonic band gap, Q± gives the upper limit for the total g-factor G f 
the cavity mode. 

Given a PhC cavity or waveguide, we can compute the near-field using Finite Difference 
Time Domain (FDTD) analysis. As described in Reference [3], the near-field pattern at a surface 
S above the PhC slab contains the complete information about the out-of-plane radiation losses 
of the mode, and thus about Q± (Fig. 1). 




Fig. 1. Estimating the radiated power and from the known near field at the surface S 
The total time-averaged power radiated into the half-space above the surface S is: 
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where K(0,<j>) is the radiated power per unit solid angle. In the appendix, we derive a very 
simple form for K in terms of 2D Fourier Transforms (FTs) of H z and E z at the surface S, after 
expressing the angles 9, in terms of k x and k y : 
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Here, T] = yj A is the mode wavelength in air, k = 2n/X, and k\\ = (k x ,k y ) = 

k(sin cos 0, sin 6 sin0) and k z = kcos(O) denote the in-plane and out-of -plane ^-components, 
respectively. In Cartesian coordinates, the radiated power (5) can thus be re-written as the inte- 
gral over the light cone, ku < k. Substituting (6) into (5) gives 
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This is the simplified expression we were seeking. It gives the out-of-plane radiation loss as 
the light cone integral of the simple radiation term (6), evaluated above the PhC slab. Substi- 
tuting Eq. (2) and Eq. (7) into Eq. (1) thus yields a straightforward calculation of the Q for a 
given mode. In the following sections, when considering the qualitative behavior of (7), we will 
restrict ourselves to TE-like modes, described at the slab center by the triad (E x ,E y ,H z ), that 
have H z even in at least one dimension x or y. For such modes, the term \FT2(H Z )\ 2 in (7) just 
above the slab is dominant, and \FT2,(E Z )\ 2 can be neglected in predicting the general trend of 

Q. 

The figure of merit for cavity design depends on the application: for spontaneous emission 
rate enhancement through the Purcell effect, one desires maximal Q/V; for nonlinear optical 
effects Q 2 /V; while for the strong coupling regime of cavity QED, maximizing ratios g/x ~ 
Qj\fV and g/y^ 1/VV is important. In these expressions, V is the cavity mode volume: 
V = (J e(r)\E(r)\ 2 d^r)/ max(£(r)|£(r)| 2 ), g is the emitter-cavity field coupling, and K and 7 
are the cavity field and emitter dipole decay rates, respectively [2]. 

Thus, for a given mode pattern, we have derived a simple set of equations that allow easy 
evaluations of the relevant figures of merit. In the next section, we address the problem of 
finding the field pattern that optimizes these figures of merit and later derive the necessary PhC 
to support it. 

3. Optimum fc-space field distribution of the cavity mode field 

The photonic crystals considered here are made by periodically modulating the refractive index 
of a thin semiconductor slab (see Fig. 2). The periodic modulation introduces an energy band- 
structure for light in two dimensions f = (x,y). The forbidden energy bands are the source of 
DBR confinement. 

The periodic refractive index e(r) = e(r + R) can be expanded in a Fourier sum over spatial 
frequency components in the periodic plane: 

e(?) = £%e^ (8) 

G 

Here G are the reciprocal lattice vectors in the (k x ,k y ) plane and are defined by G- R = 2nm 
for integer m. The real and reciprocal lattice vectors for the square and hexagonal lattices with 
periodicity a are: 



Square Lattice: 



Hexagonal Lattice: 



R m j = max + jay (9) 
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R mj = ma Vja (10) 



Fig. 2. Left: TE-like mode band diagram for the square lattice photonic crystal. The inset 
on the right shows the reciprocal lattice (orange circles), the first Brillouin zone (blue) and 
the irreducible Brillouin zone (green) with the high symmetry points. The inset on the left 
shows the square lattice PhC slab and relevant parameters: periodicity a, hole radius r, and 
slab thickness d. The parameters used in the simulation were r/a = 0.4, d/a = 0.55, and 
n = 3.6. Right: Band diagram for TE-like modes of the hexagonal lattice PhC. The inset 
on the right shows the reciprocal lattice (orange), the Brillouin zone (blue), and the first 
irreducible Brillouin zone (green) for the hexagonal lattice photonic crystal. The inset on 
the left shows the hexagonal lattice PhC slab. The parameters used for the simualtion were 
r/a = 0.3,d/a = 0.65, n = 3.6. A discretization of 20 points per period a was used for both 
diagrams. 
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where m, j, q and I are integers. The electromagnetic field corresponding to a particular wave 
vector k inside such a periodic medium can be expressed as [12]: 

G 



By introducing linear defects into a PhC lattice, waveguides can be formed. Waveguide 
modes have the discrete translational symmetry of the particular direction in the PhC lattice 
and can therefore be expanded in Bloch states Er, In Fig. 3, we plot the dispersion of a waveg- 
uide in the FJ direction of a hexagonal lattice PhC. 

Cavity modes can then be formed by closing a portion of a waveguide, i.e., by introducing 
mirrors to confine a portion of the waveguide mode. In case of perfect mirrors, the standing 
wave is described by H = a^H^ + a-^H^. (Here we focus on TE-like PhC modes, and 
discuss primarily H z , although similar relations can be written for all other field components). 
Imperfect mirrors introduce a phase shift upon reflection; moreover, the reduction of the dis- 
tance between the mirrors (shortening of the cavity) broadens the distribution of k vectors in 
the mode to some width Ak. The H : field component of the cavity (i.e., a closed waveguide) 
mode can then be approximated as: 

* +M/2 _ _ _ 

ff*(*oO~I I {A l8 e ik UA_ w e- ik - 7 )e iG ' (12) 

G k -M/2 

A similar expansion of H z can be made at the surface S directly above the PhC slab (see 
Fig. 1), which is relevant for calculation of radiation losses. The Fourier transform of the above 
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Fig. 3. Left: Band diagram for hexagonal waveguide in TJ direction, with r/a = 0.3, d/a = 
0.65, n = 3.6. The bandgap (wedged between the gray regions) contains three modes. Mode 
B oa can be pulled inside the bandgap by additional neighbor hole tuning. Right: B- of 
confined modes of hexagonal waveguide. The modes are indexed by the B-field's even ("e") 
or odd ("o") parities in the x and y directions, respectively. The confined cavity modes B „, 
B ee , and B eo required additional structure perturbations for shifting into the bandgap. This 
was done by changing the diameters of neighboring holes. 



equation gives the A:-space distribution of the cavity mode, with coefficients Aj, g and A j g. The 

distribution peaks are positioned at ±k„ ± G, with widths directly proportional to Ak. The k- 
space distribution is mapped to other points in Fourier space by the reciprocal lattice vectors G. 
To reduce radiative losses, the mapping of components into the light cone should be minimized 
[3]. Therefore, the center of the mode distribution ko should be positioned at the edge of the 
first Brillouin zone, which is the region in A:-space that cannot be mapped into the light cone 
by any reciprocal lattice vector G (see Fig. 2). For example, this region corresponds to ko = 
±jx%k x ± jyjky for the square lattice, where j x ,j y e 0,1. Clearly, \j x \ = \j y \ = 1 is a better 
choice for since it defines the edge point of the 1st Brillouin zone which is farthest from the 
light cone. Thus, modes centered at this point and £-space broadened due to confinement, will 
radiate the least. Similarly, the optimum ko for the cavities resonating in the TJ direction of the 
hexagonal lattice is fco = ±f k x (as it is for the cavity from Ref. [4], and for the FX direction is 
ko = ±^k y (as it is for the cavity from Refs. [2] and [3]). 

Assuming that the optimum choice of ko at the edge of the first Brillouin zone has been 
made, the summation over G can be neglected in Eq. (12), because it only gives additional 
Fourier components which are even further away from the light cone and do not contribute to 
the calculation of the radiation losses. In that case, we can express the field distribution as: 

H z (x,y) = J jdk x dk y A{k x ,k y )e il7 , (13) 

where A(k x ,k y ) is the Fourier space envelope of the mode, which is some function centered at 
ko = (±kox,±ko y ) with the full-width half-maximum (FWHM) determined by Ak = (Ak x ,Ak y ) 
in the k x and k y directions respectively. Eq. (13) implies that H z (x,y) and A(k x , k y ) are related by 
2D Fourier transforms. For example, if A(k x ,k y ) can be approximated by a Gaussian centered 



at ko = (ko x ,koy) and with the FWHM of (Ak x ,Ak y ), the real space field distribution H z (x,y) 
is a function periodic in the x and y directions with the spatial frequencies of ko x and koy, 
respectively, and modulated by Gaussian envelope with the widths Ax ~ 1 /Ak x and Ay ~ 1 /Ak y . 
Therefore, the properties of the Fourier transforms imply that the extent of the mode in the 
Fourier space Ak is inversely proportional to the mode extent in real space (i.e., the cavity 
length), making the problem of Q maximization even more challenging when V needs to be 
simultaneously minimized. This has already been attempted in the past for a dipole cavity 
[2, 3] and a linear defect [4, 9, 11], by using extensive parameters space search. In the following 
sections, we will design high Q cavities by completely eliminating the need for parameter space 
searches and iterative trial and error approaches. 

There are two main applications of Eq. (7). First, this formulation of the cavity Q factor 
allows us to investigate the theoretical limits of this parameter and its relation to the mode 
volume of the cavity. Second, it allows us to quantify the effect of our perturbation on the 
optimization of Q using only one or two layers of the computational field and almost negligible 
computational time compared to standard numerical methods. We applied Eq. (7) to cavities 
obtained from an iterative parameter space search. These cavities were previously studied in 
[3] and [4]. The results for Q using Eq. (7) at S, as well as full first-principle FDTD simulations 
are shown in Fig. 4. A good match is observed. Therefore, our expression (7) is a valid measure 
of the radiative properties of the cavity and can be used to theoretically approach the design 
problem; we can also use this form to speed up the optimization of the cavity parameters. The 
discrepancy between Eq. (7) and FDTD is primarily due to discretization errors. 
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Fig. 4. Comparison of Q factors derived from Eq. (7) (squares) to those calculated with 
FDTD (circles). Left: cavity made by removing three holes along the TJ direction confining 
the B oe mode. The Q factor is tuned by shifting the holes closest to the defect as shown by 
the red arrow. The x-axis gives the shift as a fraction of the periodicity a. Right: the X 
dipole cavity described in [3]. The Q factor is tuned by stretching the center line of holes 
in the FX direction, as shown by the arrow. The x-axis gives the dislocation in terms of the 
periodicity a. 



4. Inverse problem approach to designing PhC cavities 

In the inverse approach, we begin with a desired in-plane Fourier decomposition of the resonant 
mode, FT2(H(r)), chosen again to minimize radiation losses given by Eq. (7). The difficulty 
lies with designing a structure that supports the field which is approximately equal to the target 
field, H(r). 

In this section, we first estimate the general behavior of Q/V for structures of varying mode 
volume. Then we present two approaches for analytically estimating the PhC structure e(r) 
from the desired £-space distribution FT2(H(r)). As mentioned in Sec. 2, we restrict the analysis 



to TE-like modes B eo ,B oe , and B ee (Fig. 3(b)) for which we can approximate the trend of the 
radiation (7) by considering only H z at the surface S just above the PhC slab. Moreover, to make 
a rough estimate of the cavity dielectric constant distribution from the desired H z field on S, we 
approximate that H z at S is close to H z at the slab center. 

4.1. General trend of Q/V 

The simplification described above allows us to study the general behavior of Q/V for a cavity 
with varying mode volume. Here, we assume that a structure has been found to support the 
desired field H z . 

We again start from the expression for radiated power, Eq. (7), and calculate^ using Eq. (1). 
All that is required of the cavity field is that its FT at the surface S above the slab be distributed 
around the four points = ±7z/a,k y Q = to minimize the components inside the light 

cone. As an example, we choose a field with a Gaussian envelope in Fourier space, as described 
in Sec. 3. For now, let us consider mode symmetry B oe . The Fourier Transform of the H z field 
is then given by 

FT 2 (H Z )= £ sign(k M )ex V {-(k x -k M ){a x /V2) 2 -(k y -kyo){a y /V2) 2 ), (14) 

where O x and <7 V denote the modal width in real space. The mode and its FT are shown in Fig. 5 
(d-e). We use Eq.7 without the E z terms to estimate the trend in Q, as described above. As the 
mode volume grows, the radiation inside the light cone shrinks exponentially. This results in 
an exponential increase in Q. This relationship is shown in Fig. 5(g) for field B oe at frequency 
a/A = 0.248. At the same time, the mode volume grows linearly with a x . The growth of Q/V 
is therefore dominantly exponental , and we can find the optimal Q for a particular choice of 
mode volume (i.e. a x ) of the Gaussian mode cavity. 

According to Fig. 5(g), very large Qs can be reached with large mode volumes and there does 
not seem to be an upper bound on Q±. As the mode volume of the Gaussian cavity increases, 
the radiative Fourier components vanish exponentially, but are never zero. A complete lack of 
Fourier components in the light cone should result in the highest possible Q. As an example of 
such a field, we propose a mode with a sine envelope in x and a Gaussian one in y. The FT of 
this mode in Fig. 5(b) is described by 

FT 2 (H Z )= £ exp(-{k y -k }i) ) 2 (a y /V2) 2 )Rect(k x -k x0 ,Ak x ), (15) 

where Rect(k x ,Ak x ) is a rectangular function of width Ak x and centered at k x . 

The Fourier-transform implies that the cavity mode is described by a sine function in x whose 
width is inversely proportional to the width of the rectangle Rect(k x ,Ak x ). To our knowledge, 
this target field has not been previously considered in PhC cavity design. This field is shown in 
Fig.5(a-c). Though it has no out-of-plane losses, this field drops off as ~ and therefore requires 
a larger structure than the Gaussian field for confinement. 

Over the past years, many new designs with ever-higher theoretical quality factors have been 
suggested [11]. In light of our result that Q/V increases exponentially with mode size, these 
large Qs are not surprising. On the other hand, direct measurements on fabricated PhC cavities 
indicate that Qs are bounded to currently ~ 10 4 by material absorption and surface roughness 
[13, 14], while indirect measurements using waveguide coupling indicate values up to 6 • 10 5 [9, 
11]. 

It is interesting to note that Eq. (6) also allows one to calculate the field required to radiate 
with a desired radiation distribution. For example, many applications require radiation with a 
strong vertical component; waveguide modes with even H z can be confined for this purpose so 




Fig. 5. Idealized cavity modes at the surface S above the PhC slab; all with mode volume 
~ (A/n) 3 . (a-c) Mode with sine and Gaussian envelopes in x and y, respectively: H z (x,y), 
FTl(H z ), and K(k x ,k y ) inside the light cone; (d-g) Mode B oe with Gaussian envelopes in 
the x and y directions : H z (x,y), FTi(H z ), K(k x ,k y ), and (2x(<T x /a) as well as Q±/V . 
2x was calculated using Eq. (7) and E z was neglected.; (h-i) Mode B ee with Gaussian 
envelopes in x and y can be confined to radiate preferentially upward. 

that K(k Xl ky) dominates losses at the origin in k— space, as shown for instance in Fig. 5(h-i) for 
the confined mode pattern B ee of Fig. 3. 

4.2. Estimating Photonic Crystal Design from k-space Field Distribution 

Now we introduce two analytical ways of estimating the dielectric structure e(r) that supports a 
cavity field that is approximately equal to the desired field H c . These methods directly calculate 
the dielectric profile from the desired field distribution, without any dynamic tuning of PhC 
parameters, and are thus computationally fast. We focus on TE-like modes, since they see a 
large bandgap and exhibit electric-field maxima at the slab center. For TE-like modes, H c = 
H c z at the center of the slab, and H c m H c z at the surface. First, we relate H c to one of the 
allowed waveguide fields H w . The fields H c and H w at the center of the PhC slab (z = 0) are 
solutions to the homogeneous wave equation with the corresponding refractive indices e c and 
e w , respectively. 

d 2 H 1 
-Mo^ £ = « c 2 Ato4 = Vx-VxH ( . (16) 



-jUo^- = (o; v plqH w = V X —V X # w 



(17) 



Here O c and co w are the frequencies of the cavity and waveguide fields. We expand the cavity 
mode into waveguide Bloch modes: 



H c = £c,«,(r)e^-V) 
k 



(18) 



where Uk is the periodic part of the Bloch wave. Assuming that the cavity field is composed 
of the waveguide modes with k f=a ko, we can approximate u^r) f=a u^if), which leads to the 
slowly varying envelope approximation: 



H c ~ Uk (r)e 



i(k -r-(O„t) 



i((k-k )-r-(oi k -m lr )t) 



(19) 



where the waveguide mode H w = M< :o (r)e i ^ :o ' r_ ' "'' and the cavity field envelope H e = 

^ CA:e i'(^-^o)-r-(%-(«.v». 

The cavity and waveguide fields FT-distributions are concentrated at the edge of the Brillouin 
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and a <C 1 (i.e., the band is nearly flat). Differentiating (19) 
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co 2 + a\k~k Q \' 



(0 2 ,H w H e + aH w V 2 H e 



(20) 



Thus, for a <C 1 and finite V 2 H e , C0 C ~ (Ow, i.e., the cavity resonance is very close to the 
frequency of the dominant waveguide mode. The condition acl, also implies that Cfy. w CGw, 
i.e. H e «£|C*e'' ( *"*°K 

4.2.1. Estimating Photonic Crystal Design from &-space field Distribution: Approach 1 

Let us express the cavity dielectric constant e c as e c = £„ ,E e , where E e is the unknown envelope. 
H c is a solution of Eq. (16) and each waveguide mode ui c (r)e'( k ' r ~ akt ^ satisfies Eq. (17). From 
previous arguments, for k within the range corresponding to the cavity mode, cq c ~0) w k> a>£. 
Thus, a linear superposition of waveguide modes £j c^u^^e'^ : ' r ~°V) = H e H w = H c also satis- 
fies Eq. (17), i.e. the cavity mode is also a solution of Eq. (17) for a slowly varying envelope. 
We assume that the mode is TE-like, so that the H field only has a z component at the cen- 
ter of the slab. Then inserting H c from (19) into Eq.(16) and Eq. (17) and subtracting the two 
equations with E c = £ K £ e , yields a partial differential equation for £ c (x,y): 



id- 



■ lMs{col-tol)H c ttO 



(21) 



For this approach, we consider waveguide modes with B oe symmetry in Fig. 3(b). These 
modes are even in y, so the partial derivatives in y in (21) vanish at y = 0. The resulting sim- 
plified first-order differential equation in 1 /e e can then be solved directly near y = 0, and the 
solution is -?-(-! — l)d x H c w C , where C is an arbitrary constant of integration. In our analysis 



e w corresponds to removing holes along one line (x-axis) in the PhC lattice. The cavity is cre- 
ated by introducing holes into this waveguide, which means that J — 1 > 0. The solution holds 
when we take the absolute value of both its sides, and for C > 0, this leads to the following 
result for the cavity dielectric constant near y = 0: 



e c {x) 



\d x H c \ 



C+j-\d x H c 



(22) 



C is a positive constant of integration, and H c = H w H e , where H w is the known waveguide 
field and H e is the desired field envelope. C can be chosen by fixing the value of e c at some 
x, leading to a particular solution for e c . In our cavity designs we chose C such that the value 
of E c is close to £„. at the cavity center. To implement this design in a practical structure, we 
need to approximate this continuous e c by means of a binary function with low and high-index 
materials £/ and £/,, respectively. We do this by finding, in every period j, the air hole radius rj 
that gives the same field-weighted averaged index on the x-axis: 



ja+a/2 
ja-a/2 



(e h + (ei-e h )Rect{ja,rj)) 



E c 



2 rja+a/2 
dx = E c (x) 
J ja-a/2 



dx, 



(23) 



where E c is estimated from a linear superposition of waveguide modes as £ c o<Vx H c . We 
assume that the holes are centered at the positions of the unperturbed hexagonal lattice PhC 
holes. 

The radii r ; thus give the required index profile along the x symmetry axis. The exact shape 
of the holes in 3D is secondary - we choose cylindrical holes for convenience. Furthermore, we 
are free to preserve the original hexagonal crystal structure of the PhC far away from the cavity 
where the field is vanishing. 

To illustrate the power of this inverse approach, we now design PhC cavities that support the 
Gaussian and sinc-type modes of Eq. (14), (15). In each case, we start with the waveguide field 
B oe of Fig. 3(b) confined in a line-defect of a hexagonal PhC. The calculated dielectric structures 
and FDTD simulated fields inside them are shown in Fig. 6. The FT fields on S also show a 
close match and very little power radiated inside the light cone (Fig.6(c,f)). This results in very 
large Q values, estimated from Q± to limit computational constraints. These estimates were 
done in two ways, using first principles FDTD simulations [2], and direct integration of lossy 
components by Eq. (7). The results are listed in Table 1 and show an improvement of roughly 
three orders of magnitude over the unmodified structure of Fig. 4 with as small increase in 
mode volume. Furthermore, a fit of the resulting field pattern to a Gaussian envelope multiplied 
by a Sine, yielded a value of a x / a m 1 .6, which, according to the plot in Fig. 5 g., puts us at the 
attainable limit of Q± at this mode volume. 

In our FDTD simulations, we verified that Q± correctly estimates Q by noting that Q± 
did not change appreciably as the number of PhC periods in the x— and y— directions, N x 
and N y , was increased: for the Gaussian-type (sinc-type) mode, increasing the simulation size 



from N x = l3,N y = 13 (N x = 21,N y 



9) PhC periods to N x 

i3 



25, N Y 



13 (N x 

"(3 



13) 



10 6 ) to 



33, N y 

changed quality factors from g|| = 22 • 10 J ,g± = 1.4 • 10° (Q\\ = 17 • 10 J ,£>± = 4.2 " 
Q\\ = 180 10 3 ,e ± = 1.48 10 6 (gy =260- 10 3 ,g ± =4.0- 10 6 ). (The number of PhC periods in 
the x-direction in which the holes are modulated to introduce a cavity is 9 and 29 for Gaussian 
and sine cavity, respectively, while both cavities consist of only one line of defect holes in the 
y-direction.) Thus, with enough periods, the quality factors would be limited to Q±, as summa- 
rized in the table. In the calculation of Q, the vertically emitted power (Pn) was estimated from 
the fields a distance ~ 0.25 • X above the PhC surface. Note that the frequencies a/X closely 
match those of the original waveguide field B oe (a/ X cav =0.251), validating the assumption in 
the derivation. 



Table 1 . Q values of structures derived with inverse-approach 1 







Qcav (freq. filter) 


Qcav (Eq. (7)) 


Vtnode ( n ) 


Gaussian 


0.248 


1.4- 10 6 


1.6 10 b 


0.85 


Sine 


0.247 


4.2 • 10 6 


4.3 -10 6 


1.43 


Unmodified 3 -hole defect 


0.251 


6.6- 10 3 


6.4 -10 3 


0.63 
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Fig. 6. FDTD simulations for the derived Gaussian cavity (a-c) and the derived sine 
cavity (d-f). Gaussian: (a) B z ; (b) \E\; (c) FT pattern of B z taken above the PhC 
slab (blue) and target pattern (red). Sine: (d) B z ; (e) |£|; (f) FT pattern of B z taken 
above the PhC slab (blue) and target pattern (red) (The target FT for the sine cavity 
appears jagged due to sampling, since the function was expressed with the reso- 
lution of the simulations). The cavities were simulated with a discretization of 20 
points per period a, PhC slab hole radius r = 0.3a, slab thickness of 0.6a and refrac- 
tive index 3.6. Starting at the center, the defect hole radii in units of periodicity a are: 
(0,0,0.025,0.05,0.075,0.1,0.075,0.075,0.1,0.125,0.125,0.125,0.1,0.125,0.15,0.3,0.3) 
for the sine cavity, and (0.025,0.025,0.05,0.1,0.225) for the Gaussian cavity. 
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4.2.2. Estimating Photonic Crystal Design from &-space Field Distribution: Approach 2 

We now derive a closed-form expression for e c (x,y) that is valid in the whole PhC plane (instead 
of the center line only). Again, begin with the cavity field H (r) = zH c consisting of the product 
of the waveguide field and a slowly varying envelope, H c = H w H e , and treat the cavity dielectric 
constant as: -J- = — ■ — h 4~- In the PhC plane, Eq. (16, 17) for a TE-like mode can be rewritten 

as 

-0) c 2 AioH, = V-(ivH c ) (24) 
-fl^Moffw = VH W ) (25) 

Multiplying the last equation by H e , subtracting from the first, and recalling that co c ~ C0 W yields 

collM)H e H w -co^oHc = /ioff c (a>jJ-fl£)M() (26) 
= V-(-VH c )-H e V-(—VH w ) 

« V-(— VH C ) (27) 

£pert 



where the last line results after some algebra and dropping spatial derivatives of the slowly 
varying envelope H e . This relation is a quasilinear partial differential equation in 1 /B pert . With 
boundary conditions that can be estimated from the original waveguide field, this equation can 
in principle be solved for e c (e.g., [15]). 

Alternatively, one can find a formal solution for e pe rt by assuming a vector function t, (?) 
chosen to satisfy the boundary conditions, so that 

—VH C = V x | (28) 



tpert 

or 

1 V x | • VH* 



(29) 

pert | v ll c\ 

This gives is the formal solution of the full dielectric constant e c — (SpJ^ + e" 1 ) _1 in the plane 



solution of the full diele 

of the photonic crystal. 



5. Conclusions 

We have described a simple recipe for designing two-dimensional photonic crystal cavities. Al- 
though the approach is general, we have demonstrated its utility on the design of cavities with 
very large Q > 10 6 and near-minimal mode volume ~ (A /«) 3 . These values follow our theoret- 
ically estimated value of Q±/V for the cavity with the Gaussian field envelope, which means 
that we were able to find the maximal Q for the given mode volume V under our assumptions. 
Our approach is analytical, and the results are obtained within a single computational step. We 
first derive a simple expression of the modal out-of-plane radiative loss and demonstrate its 
utility by the straightforward calculation of Q factors on several cavity designs. Based on this 
radiation expression, the recipe begins with choosing the FT mode pattern that gives the desired 
radiation losses. For high-Q cavities with minimal radiative loss inside the light cone, we show 
that the transform of the mode should be centered at the extremes of the Brillouin Zone, as far 
removed from the light cone as possible. Next we proved that for a cavity mode with a Gaussian 
envelope, Q/V grows exponentially with mode volume V, while the cavity with the sine field 
envelope should lead to even higher Q's by completely eliminating the Fourier components in 
the light cone. Finally, we derived approximate solutions to the inverse problem of designing a 
cavity that supports a desired cavity mode. This approach yields very simple design guides that 
lead to very large Q/V. Since it eliminates the need for lengthy trial-and-error optimization, 
our recipe enables rapid and efficient design of a wide range of PhC cavities. 
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A. Derivation of Cavity Radiative Loss 

The radiated power per unit solid angle K(9,(j>) can be expressed in terms of the radiation 
vectors N and L in spherical polar coordinates (r, , ) : 



*(0,0) = 



1 

8A 2 



L<j> 



(30) 



where n = . /■f 2 . The radiation vectors in spherical polar coordinates can be expressed from 

y *~o 

their components in Cartesian coordinates: 



Ne = (A^cos^ +Af y sin0)cos 9 
Afy = -N x sin (j> +N y cos (j>, 



(31) 



and similarly for Lg and L§. As described in Reference [3], the radiation vectors in Cartesian 
coordinates are proportional to 2D Fourier transforms of the parallel (x and y) field components 
at the surface S (Fig. 1): 



N x 

Ny 

Lx 

Ly 

h 
k 7 



-FT 2 (H y ) 
FT 2 (H X ) 
FT 2 (E y ) 
-FT 2 (E X ) 



x v 

k( — , — ) = k sin 9 (cos ((>x + sin (j>y) 



(32) 



A: cos 9, 



where k = 2n/X, X is the mode wavelength in air, and kn = A sin 9. 
Here the 2D Fourier Transform of the function f(x,y) is 



FT 2 (f(x,y)) = J J d x dyf{x,y)e fk W^ (33) 
= J jd.jlJ>x.y)c n ^ k -< (34) 

Substitution of expressions (32) and (33) into (30) now yields an expression for the radiated 
power (5) in terms of the FTs of the four scalars H x ,H y ,E x , and E y . This expression is in general 
difficult to track analytically. We will now use Maxwell's relations to express (30) in terms of 
only two scalars, H z and E z . 

Noting that for a bounded function g, FT 2 (^) = —ik x FT 2 (g) (similarly for ^), we can 
re-write Ng as 

N e = -j^(-k x FT 2 (Hy)+k y FT 2 (H x )) (35) 

. k z pT dHy dH x ^ 
k\\k " dx dy 

= —j—FT 2 (E z ) 

where the last step follows from Maxwell's Eq. V X H = £i>~jf = iC0£ o E. Similarly, we find 
Lg = ^FT 2 (H Z ). From V • H = and V • E = at the surface S, it also follows that N$ = 



T-FT 2 (^) and = j-FT 2 {^). Substituting these expressions into Eq. (30) yields 



K(k x ,k y ) 



U 2 k 2 



k z FT 2 (E z ) + iFT 2 (^-) 



k z FT 2 {H z )+iFT 2 (-^) 
dz 



(36) 



Furthermore, E Z ,H Z °< e\p(ik z z) for propagating waves inside the light cone (which are the 
only ones that determine P), implying that FT 2 {^t) = —ik z FT 2 (E z ) and similarly for H z . This 
allows further simplification of the previous expression to 



K(k x ,ky) 



x\k\ 



'FT 2 (E Z )\ 2 + \FT 2 (H Z )\ 2 



(37) 



Substituting this result back into the expression for total radiated power (5) gives the required 
result (7): 



J J d9d((>sm{e)K(e,<l>) 



(38) 



o o 



dk x dk y —K(k Xl k Y ) \J{k x , k y ) 

k«<k ' k 



1 



2X 2 kJk n <k k 



2 



\\FT 2 (E Z )\ 2 + \FT 2 (H Z )\ 2 



Above, J(k x ,ky) is the Jacobian resulting from the change of coordinates from {9,(j)) to 
(k x ,ky). 



